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Abstract 

Massive parallelisation has lead to a dramatic increase in available 
. £h computational power. However, data transfer speeds have failed to 
Ay keep pace and are the major limiting factor in the development of 
J3 exascale computing. New algorithms must be developed which minimise 
the transfer of data. Patch dynamics is a computational macroscale 
modelling scheme which provides a coarse macroscale solution of a 
problem defined on a fine microscale by dividing the domain into 
many nonoverlapping, coupled patches. Patch dynamics is readily 
adaptable to massive parallelisation as each processor can evaluate 
the dynamics on one, or a few, patches. However, patch coupling 
conditions interpolate across the unevaluated parts of the domain 
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between patches, and are typically reevaluated at every microscale 
time step, thus requiring almost continuous data transfer. We propose 
a modified patch dynamics scheme which minimises data transfer by 
only reevaluating the patch coupling conditions at ‘mesoscale’ time 
scales which are significantly larger than the microscale time of the 
microscale problem. We analyse the error arising from patch dynamics 
with mesoscale temporal coupling as a function of the mesoscale time 
interval, patch size, and ratio between the microscale and macroscale. 
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1 Introduction 

Mathematical equations describing a phenomenon (e.g., fluid flow, chemotaxis, 
mechanics) are typically written at the scale at which we want the informa¬ 
tion (e.g., macroscopic velocity fields, bacterial concentrations, macroscopic 
deformations). Increasingly, the scale at which the physics are understood 
(molecular, cellular, agent-based) is much finer than the macroscopic, human, 
systems scale at which we want information and a variety of modelling tech¬ 
niques may be applied to reinterpret the microscale problem at the desired 
macroscale [3, 12, 15, 8, 9, e.g.]. Unfortunately, for many multiscale and 
multiphysics problems good macroscale closures do not exist—instead we 
simulate and observe very detailed models of great complexity at great cost 
[26, 25, 20, 24, e.g.]. In this scenario, we aim to develop efficient computa¬ 
tional procedures to wrap around whatever microscopic level computer model 
a scientist chooses for any given system [18, 17, e.g.]—be it anything from a 
Monte-Carlo description of a chemical reaction to an individual/agent-based 
model in ecology or epidemiology. 

The methodology is to evaluate automatically (‘on demand’), directly from 
the micoscale model, the macroscopic modelling closures for the emergent- 
dynamics which all too often are not available explicitly [7, e.g.]. This is 
an ‘equation-free’ method in the sense that it makes no attempt to derive a 
macroscale equation, in contrast to, for example, homogenization [27]. 

The aim of this article is to develop and support the patch dynamics 
scheme [12, 16, 22, 19, 32, for reviews]. By only computing on a small 
fraction of the space-time domain, see the schematic Figure 1, this scheme 
empowers large scale simulation and prediction. But special challenges arise 
in the largest simulations on exa/pet-a-scale computers. Designed for exa/ 
pet-a-scale computing we propose new infrequent couplings between these 
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Figure 1: Equation-free methods simulate only on small patches in space-time: 
patches are placed at macroscale time steps At and spatial macroscale steps H. 
The given microscale dynamics are only simulated within each patch: coupling 
conditions interpolate across un-simulated space; and projective integration 
steps across time. 

microscale simulations on microscale patches across un-simulated space. We 
establish new results on efficiency, accuracy, and consistency for the emergent 
macroscale simulation. Section 2 discusses the mathematical details of patch 
dynamics with patches defined in space only (also referred to as the gap-tooth 
method [30]) and uses the example of a one dimensional diffusion problem on 
a discrete lattice. 

The scheme is readily adaptable to higher dimensions and more complex non¬ 
linear systems. Section 6 develops numerical simulations of a two dimensional 
complex Ginzburg-Landau pde. Simulations, such as that shown in Figure 2, 
qualitatively confirm the accuracy of our proposed infrequent coupling scheme 
for this two dimensional nonlinear system. 

As computational power approaches the exascale, roughly 10 18 FLOPS (float¬ 
ing point operations per second), it is tempting to think that soon many 
multiscale problems will be solved numerically by computing a microscale 
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u(x,y,t) at t = 00.0400 


u(x, y, t) at t = 00.4000 



Figure 2: Real parts of Ginzburg-Landau microscale fields u^+^N^+iyN * n 
patches with n = 6 at time t = 0.04 and t = 0.4 with: (top) continuous time 
coupling; and (bottom) infrequent mesoscale coupling 6t = 0.2. At this scale 
there is little to distinguish the continuous time coupling solution and the 
mesoscale coupling solution. 
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simulation across the entire domain. However, constraints on high perfor¬ 
mance computing make such a task effectively impossible for all but simple 
scenarios [10]. Improvements in high performance computing are gained 
through massive parallelisation via increases in the number of processors, but 
success is forecast to be limited [21, 35, 2, 37, 11]. First, one limitation is 
that memory storage is growing at a tenth of the pace of processing power [2], 
so there is huge processing power but little space to store the resulting data. 
Schemes for exascale computing must use only as much data as required; that 
is, we need as sparse as possible a resolution over the macroscale, such as that 
offered by the patch scheme. Secondly, relative to computational speeds, the 
slow speed of data transfer between processors, cache and memory prevents 
processors from operating effectively unless the computational scheme limits 
communication, as we propose and analyse here for the patch scheme. Thirdly, 
with millions of processors, hardware failure will be common somewhere and 
we discuss possibilities for fault tolerance in the patch scheme. These limita¬ 
tions are not expected to be overcome through improvements in hardware, so 
it is mainly through the development of new algorithms that engineers and 
scientists will exploit the benefits of massive parallelisation [36, 2, 37, 11]. 

The patch dynamics approach does not invoke a macroscale equation and 
requires no prior analysis of the spatial-temporal domain [19, e.g.]. This 
removes significant data storage constraints while increasing the flexibility, 
allowing ‘on-the-fly’ modifications. The discretisation of the domain into 
patches makes patch dynamics readily adaptable to massive parallelisation: 
for example, a domain decomposition where each processor simulates the 
dynamics on a few patches. Here, for simplicity, we assume only one patch per 
processor. However, extant implementations of patch dynamics require that 
coupling conditions are calculated and communicated between each patch at 
each microscale time step [16, 28, 29, 30, 31], thus requiring substantial and 
effectively prohibitive data transfer between processors. Section 2 proposes a 
new modification to the patch scheme to reduce data transfers for exascale 
computing by limiting the times at which the inter-patch coupling conditions 
are updated. As illustrated by the red and blue lines in Figure 3, coupling 
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condition data required from other patches (or processors) is updated only at 
mesoscale time steps St which are significantly larger than the microscale time 
steps of the simulator but smaller than the macroscale time of interest At. 
However, as indicated by the brown line in Figure 3, the data required for 
one patch’s coupling conditions which is dependent on the dynamics of that 
patch is updated at microscale time intervals since this information is readily 
available to the processor. These mesoscale coupling adjustments to the patch 
dynamics scheme should greatly increase the speed of a simulation run on a 
high performance computer. 

A significant issue for exascale computing is fault management and re¬ 
siliency [36, 37, 11]. Typically, if a computer component fails and data 
is lost, then, assuming the computer is still operational, the required calcula¬ 
tion is redone either from scratch or, if there is some fault tolerance written 
into the algorithm, from a checkpoint. If a failure causes data to be delayed 
rather than lost, then the whole computation is delayed until the data is 
successfully transferred. In either case, failures increase the time required 
to complete a calculation. Such delays are not usually major issues for a 
computer with relatively few components, but on an exascale computer with 
millions of processors and numerous other components, failure is expected to 
be a regular occurrence and restarting from a checkpoint and waiting for data 
is not viable [34, 2, 11]. Algorithms for exascale computing must be fault 
tolerant while also accounting for errors associated with failure. Section 7 
discusses how fault management may be incorporated into the proposed patch 
dynamics scheme. 

To demonstrate how to apply patch dynamics mesoscale coupling, Section 3 
solves a fundamental microscale discrete diffusion problem using standard 
patch dynamics macroscale modelling without mesoscale coupling; that is, 
with inter-patch coupling at microscale time intervals. Then, Section 4 
modifies the solution to allow for patch coupling at only mesoscale time 
intervals. Section 5 analyses the error of the solution with mesoscale coupling 
obtained in Section 4 relative to the solution obtained without mesoscale 
coupling obtained in Section 3. Therefore, this error is not the full error of the 
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Figure 3: Schematic description of the coupling conditions (5) or (6) for the 
ith patch with nearest neighbour coupling where coupling between patches 
is only reevaluated at mesoscale time steps St. Three patches are indicated 
by shaded regions centred about macroscale lattice points Xi_i , and Xi+i . 
The macroscale lattice spacing is H. As indicated by the coloured lines, 
the average dynamics on patches I ± 1 and i (2) feed into the coupling 
conditions (5) or (6) on the ith patch and this controls the dynamics on both 
edges of the ith patch. At time mSt + t, for nonnegative integer m and 
At < t < At + St, the coupling conditions of the ith patch are dependent 
on the average dynamics of the ith patch at time mSt +1 (brown) as well as 
the average dynamics of the neighbouring patches from the mesoscale time 
step mSt (blue, red). 
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modelling but requires additional consideration of the error associated with 
the standard patch dynamics scheme. The error of standard patch dynamics 
has been discussed for several different systems with a variety of microscale 
structures [28, 29, 33, 5, 31, 6]. Although we only consider the mathematical 
details for one patch, it is readily scalable to a multiple processor system, as 
shown in the numerical results of Section 6. 


2 Patch dynamics implementation 

As an initial prototype problem we consider a simple diffusion system on a 
discrete one dimensional microscale spatial lattice with lattice index j and 
microscale lattice spacing h, and time t which is measured on some microscale, 

Uj(t) =Uj+i(t) +Uj_i(t) -2uj(t), (1) 

with some given initial condition iq(0) on the microscale held. Realistic 
microscale dynamics are much more complicated than this simple microscale 
diffusion but before we can contemplate realistic dynamics we prove that 
the proposed procedure for exascale computing is sound for at least this 
foundational case of the system (1). Section 6 presents successful numerical 
simulations of the more complex two dimensional Ginzburg-Landau ODE. 

Suppose we require a macroscale simulation of ode (1) but only at discrete 
spacings H h. We construct a macroscale lattice with spacing H and 
macroscale lattice sites Xi = iH = iNH, for patch index i = 0, ±1,... and 
where the number of microscale lattice points within one macroscale step 
N = H/h. is integral. In the patch dynamics scheme, for all i we construct the 
ith patch of width 2nh, for positive integer n, centred about the macroscale 
lattice site Xi : n/N < 1 /2 ensures the patches do not overlap; in practice, 
n/N <C 1 /2 for efficient macroscale simulation. In Figure 3 the shaded areas 
schematically represent three patches at Xy-i and Xi , and in Figure 4 these 
same patches are superimposed on both the microscale and the macroscale 
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2nh 


Xt—i H X t Y X i+1 ^ 

Figure 4: The microscale lattice with microscale lattice spacing h, the 
macroscale lattice with N = 20 microscale lattice points in one macroscale 
step, and three patches with patch half-width n = 8. The relatively large 
ratio n/N = 0.4 means there is little space between the patches. 


lattices. Integer n is the patch half-width; that is, n is the number of 
microscale lattice points which fit into half a patch. Figure 4 shows patches 
with patch half-width n = 8 and n/N = 0.4. The ratio n/N is equal to the 
ratio of half the physical patch width nh. to the macroscale lattice spacing H. 

We solve ODE (1) for microscale fields Uj(t) but only for the microscale lattice 
points which lie within a patch; that is, for with microscale sub-patch 

index j = 0, ±1,..., ±(n — 1) and macroscale patch index i. Without loss of 
generality we start at t = 0, although the results presented here apply to any 
initial time which is an integer multiple of mesoscale time step 5t. We assume 
the initial and final simulation times are in the same patch, such as those 
times shown in Figure 4, so do not consider patch coupling across time. The 
patch boundary conditions of the microscale simulators are provided by patch 
coupling conditions which extrapolate across the un-simulated space between 
patches. Coupling between adjacent patches is achieved by constraining the 
average of the microscale field near both edges of each patch, termed the 
‘action regions’. The left and right action regions of one patch are shaded blue 
in Figure 5 and extend over microscale sub-patch indices j = ±n,..., (n — 2a) 
for some integer a, 0 < a < n . After solving for all microscale fields within 
the ith patch we extract the desired macroscale solution Ui(t) at Xi via 
some averaging over the microscale solutions in the middle or ‘core’ of the 
ith patch. The core of one patch is shaded brown in Figure 5 and extends 
over j = 0, ±1,..., ±a . The integer a is defined as the core half-width. Li 
et al. [23] used similar averaging techniques over action regions and core to 
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2rH = 2{n — a)h 
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Figure 5: A closeup example of the ith patch, similar to that shown in Figure 4 
with patch half-width rt = 8 . The core and action half-width is a = 2 , and 
the averaging over the core (brown) provides the macroscale held (2). The 
left and right action regions (blue) for the coupling conditions (3), (5) or (6) 
are averages over j = ±(8, 7,... ,4). The ratio r = (n — a)/N = 0.3, using 
N = 20 from Figure 4. The buffer width is n — a = 6 = rN . While it is only 
the core which contributes to the evaluation of the macroscale held Ui(t), 
extending the domain of the patch into buffer regions improves the accuracy 
of the macroscale held. 

simulate molecular dynamics in a huid. 

The helds Uj +iN in the core (i.e., j = 0, ±1,...,±a) of each patch i are 
the most useful as they dehne the macroscale solution Ut(t) which is used 
for both interpolation between patches and extrapolation across time. In 
contrast, the helds Uj+iN which lie inside a patch but outside the core (i.e., 
j = ±(a+l),..., ±n) are of no interest and are forgotten as soon as they are 
calculated. However, this region inside a patch but outside the core performs 
an important function in that it ‘sheilds’ or ‘buffers’ the helds within the 
core from errors which arise on the application of coupling conditions on 
the microscale solutions in the action regions. The left and right ‘buffers’ 
extend from the patch edge to the core edge over sub-patch indices j = 
± ( a + 1),..., ±n and have width n — a [6] . Generally a larger buffer results 
in smaller numerical error; Figure ff in Section 5 shows errors decreasing 
approximately exponentially with buffer width. Thus the buffers perform 
an important but ancillary task, improving the macroscale solution without 
directly contributing to its evaluation. Section 3 shows that n — a, rather 
than n, determines eigenmodes of the microscale solution in one patch and 
so, in addition to being the buffer width, n — a plays the role of an effective 
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patch half-width. We name n — a the reduced patch half-width. 

Patch coupling conditions and the method of deriving the macroscale solution 


may vary, depending on the problem being considered. In choosing these rules 
the form of the microscale model is important. For example, Section 5 shows 
that numerical error decreases with increasing reduced patch half-width n — a, 
so it is recommended that n — a is as large as possible. This recommendation 
is suitable for many smooth microscale models such as the simple diffusion 
problem (1) and the two dimensional Ginzburg-Landau equation discussed 
in Section 6. However, when the microscale model has some rough fine scale 
structure, such as a periodic spatial roughness, more care needs to be taken 
when choosing n — a. In the case of some periodic spatial roughness, minimal 
errors are obtained when patch dynamics adequately averages over the spatial 
structures; for example, when the period exactly divides the reduced patch 
half-width n — a [6]. Since ODE (1) has no rough microscale structure we do 
not need to consider this complication. 

The macroscale field obtained from each patch is generally some average over 
the microscale fields in the centre core of the patch. For core half-width a, 
0 < a < n , we average over the 2a + 1 microscale fields in the core of the 
ith patch to define the macroscale field of that patch as 



( 2 ) 


Figure 5 shows one patch with a core half-width a = 2. In two or more 
dimensions the core may be a rectangle or rectangular prism in the centre 
of a patch and, as in one dimension, the macroscale field of a patch is the 
average of those microscale fields within the core. Section 6 implements 
square patches with square cores. 

Each patch is coupled to its near neighbours via a control acting on the the 
boundary action regions of each patch. We define left and right action regions 
on the patch edges, both containing 2a + 1 microscale lattice points, with 
j = —n, —n + 1,..., — (u — 2a) and j = n, n — 1,..., (n — 2a) respectively. 
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With coupling at microscale times, the coupling conditions of the ith patch 
are proposed to be that the action region average [6] 

±n 


j=±n—2 q 

for some fi )= t n and cos t derived via classic Lagrange interpolation of neigh¬ 
bouring patch macroscale fields Uy-i (t), U i+b tb... and the patch macroscale 
held Ui.(t) (this interpolation has been proved to be effective for PDEs [30]). 
The right hand side of the coupling conditions (3) only contains macroscale 
fields obtained from applying averaging (2) to microscale fields within the 
core. Thus, as discussed above, microscale fields within the buffers make no 
direct contribution to the interpolation across the un-simulated space between 
patches. The action regions are the same size as the core and a is both the 
core half-width and the action half-width, as illustrated in Figure 5. 

The details of the patch coupling are contained in cos t and fi,± n which are 
functions of a parameter y, which controls the coupling strength between 
patches, and the ratio r = (n — a)/N , which compares half the physical 
reduced patch width (n — a)ft with the macroscale lattice spacing H. For 
our purposes, the details of the coupling contained in cos £ and are not 
important. Section 5 shows that errors associated with mesoscale coupling 
are not dependent on cosf and the error analysis is presented in terms of 
units of fi i= tn and its temporal derivatives. The derivation of cos £ and fi,± n 
to any order of coupling (i.e., the number of patches coupled to any one 
patch) is presented elsewhere [31, 6]. For example, for only nearest neighbour 
coupling these function have a parabolic dependence on the ratio r and linear 
dependence on coupling parameter y: cosf = (1 — r 2 y) and 

fi,±n(t) = l(2a+ l)ry[(r ± 1)U i+1 (t) + (r^l)U H (t)], (4) 

where the scaling by 2a + 1 and the subscript ±n on f i)±T1 are for later 
convenience. Since 0 < r < 1 /2 and 0 < y < 1 , for nearest neighbour 
coupling 0.75 < cost < 1 . Higher order couplings extend to next nearest 



L+iNlT) 


— i i.m 


fi -Pn ft) 


( 3 ) 


Thursday 9 th April, 2015 




2 Patch dynamics implementation 


14 


neighbours bh ±2 , and beyond, and contain terms of higher order in r and y; 
typically, cost > 0.6, even for high orders of coupling [6]. 

The physical problem of interest is at full coupling y = 1 . In this case the 
coupling condition (3) is effectively a Taylor series expansion of the macroscale 
fields about the centre of the action region j = ±rH = ±(n — a)h. (for left 
and right cases), set equal to the average of the microscale fields in the left 
and right action region [31, 6]. However, centre manifold theory provides 
full physical support for patch dynamics within some domain about y = 0 
(y = 0 is the no coupling case where patches have no influence on each 
other [29, 30, 32]). Generally we must consider the full range 0 < y < 1 when 
providing theoretical support for our patch dynamics scheme; however, as we 
here avoid formal definitions of ft^ and cos £ we do not delve into a detailed 
analysis of the limiting behaviour of coupling parameter y. 

Section 3 determines all eigenvalues and eigenvectors on an arbitrary single 
patch with patch half-width n, core half-width a, and the coupling condi¬ 
tions (3). From these we construct the microscale held solution of (1) on the 
ith patch with the coupling conditions (3). Importantly, in this solution on 
the ith patch, the form of the coupling to neighbouring patches, represented 
by f^in, and the coupling to the ith patch, represented by cos i, is arbitrary 
and thus is valid for any form of inter-patch coupling. Once the microscale 
solution on the ith patch is determined, core averaging (2) provides the 
macroscale solution Ut. 

In the standard patch dynamics scheme, the coupling conditions (3) are 
evaluated at each microscale time step as required by the microscale simulator— 
here the ODE (1). On a computer with massive parallelisation, and in 
the scenario where one processor simulates the system (1) over one patch, 
when applying coupling conditions (3) the coupling data must be transferred 
between processors each microscale time step. As discussed in Section 1, 
frequent data transfers defeat massive parallelisation. We propose to limit 
the transfer of data between processors by only communicating data between 
patches on mesoscale time-steps 6t, larger than the microscale time-steps of 
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the simulation but smaller than the macroscale times of interest At. Thus, as 
a first approximation we replace coupling conditions (3) with 


L 

j=±n—2a 


Uj +iN (mSt + t) 

2a+1 


bh(mSt + t) cos i + 


fi,±n(ra5t) 

2a+1 


( 5 ) 


where m = 0,1,...,M..0<t<St. Thus data transfers between processors 
are required much less frequently: the cost is an error in the simulation which 
we analyse in Section 5. Figure 3 illustrates this new scheme in the case of 
nearest neighbour coupling. 

More sophisticated mesoscale coupling conditions than (5) are obtained by 
approximating fi,± n (t) in coupling condition (3) with the first Q terms of its 
Taylor series expansion about the previous mesoscale time step mSt. Using 
f q to denote the qth derivative of coupling function f, this generalises the 
infrequent coupling conditions (5) to 


L 

j=±rt—2a 


Uj +iN (mSt + t) 
2a + 1 


Q-i 

Lh(m6t + t) cos l + y~ 

q=0 


(2a+1)q! ’ 


( 6 ) 


where the Taylor series expansion is (Q — 1)th order accurate. The error in 
simulations with such a Qth order coupling is also analysed in Section 5 

Section 4 modifies the microscale field solution with coupling conditions (3) 
obtained by Section 3 on the arbitrary ith patch to a solution with mesoscale 
coupling conditions (5) or (6). We use these new solutions to systematically 
explore the errors of various coupling schemes. Section 5 analyses the error 
which arises when coupling conditions (3) are replaced with the mesoscale 
coupling conditions (5) or (6) and its dependence on parameters such as the 
patch half-width n, core half-width a, ratio r, and Taylor series order of 
accuracy Q — 1. We consider both the error on the microscale lattice within 
the ith patch and the error of the macroscale solution l_h obtained from core 
averaging (2). 
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3 Microscale sub-patch dynamics 

To assess the error due to coupling at infrequent mesoscale time steps (5) 
or (6), this section establishes the microscale solution within one arbitrary 
patch with coupling at microscale time steps (3). Section 4 modifies the 
solution with microscale coupling (3) into a solution with mesoscale coupling 
(5) or (6). 

To construct the microscale solution within the ith patch we first rewrite 
equations (1) and (3) as a single matrix equation and then determine all 
eigenvalues and right and left eigenvectors. The microscale solution is a linear 
combination of terms in these eigenvalues and eigenvectors. This solution 
is valid for any patch half-width n and any equal sized action regions and 
core, 0 < a < n , with the exception of some special cases where repeated 
eigenvalues are associated with linearly dependent eigenvectors and the set 
of all eigenvectors do not form a complete basis on the patch. These special 
cases require generalised eigenvectors to provide a full microscale solution [4]; 
however, we do not consider generalised eigenvectors here as they further 
complicate the problem without providing any additional insights. 

The case when the microscale is the diffusion PDE, will be obtained as the 
limit of n —> oo with h. —> 0 and finite patch width 2nH. 


3.1 Matrix form 

In matrix form, the system of ODEs (1) within the ith patch and with coupling 
conditions (3) is 

Bu(t) = £u(t) + f(t), (7) 

where (2rt + 1) dimensional vector u = (u_ n+ iN, ... , u n+ ii\i) describes the 
field Uj + iN at every sub-patch coordinate j = —n, —n + 1,..., n — 1, n, and 
has initial condition u(0) = Uo = (u_ n+ iN(0), ... ,u n+iN (0)) . The forcing 
vector f = (ft _ n , 0,..., 0, fi iTl ) where fi,± n describes the coupling of the 
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ith patch to neighbouring patches, such as the nearest neighbour coupling (4). 
Matrices B = diag(0,1,1,..., 1,0) and C are (2n+1) x (2n+1) • Rather than 
use the usual matrix numbering of rows and columns (i.e., 1,2,..., (2n + 1)), 
we index the rows and columns of B and C as —n, —n + 1,..., n — 1,n, since 
these correspond to our patch indices j = —n, —n + 1,..., n — 1, n, and 
similarly for the indexing of components of the vectors u and f. For example, 
£_ n _ n is the element in the first row and first column of £. 

With the exception of the first and last rows, the nonzero elements of matrix C 
are £j,j_i , £j,j+i = 1 , £y = —2 for ) ^ ±n , which describes the discrete 
diffusion in ode (1). The first and last rows of C represent the patch coupling 
conditions (3) with the macroscale fields defined by the average over the core 
in equation (2). For an action and core half-width a, £± n ,j = £± n j + , 

where 

£a = f- 1 > n-2a<±j<n, £C [cost, -a<j<a, 

±n ’’ 1 0 otherwise, ] j 0 otherwise. 

The action and core half-width is only restricted by the size of the patch, 
0 < a < n. The action regions overlap the core when a > n/3, which 
means some microscale fields appear both in the average over the core (2) 
and in the averages over the action regions of the coupling conditions (3). 
Equation (8) and the subsequent microscale solution remain valid whether or 
not the action regions and core overlap. For example, when a = n/3 the core 
and action regions overlap at two microscale lattice points ±u/3 + IN and 
£±n,±n/3 = -1 + cost. 

We now solve the generalised eigenproblem and adjoint eigenproblem 

[C — A k B)v k = 0 and z k (£ — A k B) = 0 7 , (9) 

for right and left eigenvectors v k and z k , respectively, and eigenvalues A k . We 
normalise all eigenvectors such that z k Bv k / = 5 kk , for Kronecker delta S kk / 
and all k, k' = 0, 1 ,..., 2n — 2. The eigenproblem produces two set of, 
at most, (2n — 1) linearly independent right and left eigenvectors indexed 
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k = 0, 1 ,..., 2n — 2 , although the associated [2n — 1 ) eigenvalues may 
have multiplicity greater than one. A set of (2n — 1 ) linearly independent 
right or left eigenvectors forms a complete basis which spans the subspace 
of the microscale held vector u satisfying coupling conditions (3). 1 Some 
special cases are degenerate and do not provide the required (2n — 1) linearly 
independent eigenvectors: for simplicity we limit analysis to the generic case 
of (2n — 1) linearly independent eigenvectors. 


3.2 Act and sample on a lattice point 

We first consider the simplest case where there is no averaging over the action 
region and core, a = 0; instead we just act and sample the microscale held at 
end points and the mid point of a patch. The patch coupling conditions (3) 
only constrain the microscale helds Uj+iN at the patch edges, j = ±n, and 
the macroscale held value in this ith patch th is the microscale value at the 
centre of the patch, j = 0: equation (2) reduces to U t (t) = u iN (t). The only 
nonzero elements of the hrst and last rows of C are £. n> _ n = C n<n = — 1 and 
£±n,o = cos t. In this case we always obtain (2n — 1) linearly independent 
right or left eigenvectors and (2n — 1) distinct eigenvalues, and thus, for this 
case, the sets of all right or left eigenvectors forms a complete basis (we never 
need generalised eigenvectors). 

Using Matlab we evaluate the eigenvalues and eigenvectors of the matrix 
equations (9) for several n and a. From these numerical examples we de¬ 
termine the general analytic forms of the eigenvalues and eigenvectors and 
confirm by substitution into the matrix equations. The eigenvalues of the 
matrix equations (9) are 


A k = —2 [1 — cos(7rl k /2n)] , (10) 

1 The size of the basis which spans the subspace of the microscale field vector u is the 
number of elements of u (i.e., the number of microscale fields on one patch (2n + 1 )) minus 
the number of constraints on the elements of u . There are two constraints, provided by 
the two coupling conditions (3), so the size of the basis is (2u — 1). 
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for eigenvector mode indices k = 0,1,..., 2(n — 1) and corresponding wave 
number 

Jk+ 1 + (—1) k/2 (2£/7r— 1) for even k, 
k [k+1 for odd k. ( ^ 

The right and left eigenvectors for all sub-patch coordinates j, are, for even k, 

(vic)j = (rtsin£)~ 1/2 cos(jl k 7t/2ri), (12) 

(zic)j = (nsinf) _1/2 sin[£ - (-1 ) k/2 |j|l k 7r/2n], except (z k ) ±n = (zic)±(n-i)> 

(or, equivalently (z k )j = (—1 ) k/2 (risin£)~ 1/2 sin[(n — |j|)l k 7x/2n] for j ^ ±n) 
and for odd k, 

(v k )j = n _1/2 sin(jl k 7t/2n), 

(z k )j = nT 1/2 sin(jl k 7t/2n), except (z k ) ±n = (z k ) ±(n _q . (13) 

These eigenvectors are normalized such that z^Bv k / = S kk / for all k, k' = 
0,1,... ,2(n — 1). 

Figure 6 plots right and left eigenvectors associated with the four lowest 
magnitude eigenvalues, for patch half-width n = 20 and a = 0. These 
eigenvectors are only defined at discrete patch lattice points j = 0, ±1,..., ±n, 
but, for clarity, we plot curves rather than points. These plots are typical for 
any n. For odd k the right and left eigenvectors are identical sine functions 
which are odd about j = 0 and typical of what one sees in simple one 
dimensional diffusion [14]. For even k the right and left eigenvectors are even 
functions about j = 0 but, with their dependence on the parameter £, are 
not typical of eigenvectors of simple one dimensional diffusion. For all even k 
the eigenvectors’ amplitudes are no more than (nsin £)~ 1//2 ; however, in the 
centre of the patch, j = 0, the right and left eigenvectors behave differently 
with respect to £. For even k at j — 0 the right eigenvector is (nsin£) -1//2 
but the left eigenvector is (sinf/n) 1 / 2 . So, for small £ the right eigenvector is 
large at j = 0 but the left eigenvector is small, but as £ increases the right 
eigenvector at j — 0 decreases and the left eigenvector at j — 0 increases. 
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Figure 6: Right and left microscale eigenvectors associated with the four 
lowest magnitude microscale eigenvalues (k = 0,1,2,3 for blue, green, red and 
cyan, respectively) for patch half-width n = 20 , action and core half-width 
a = 0 and cost = 0.91 (so (nsin£)~ 1/2 = 0.35 and (sin £/n) 1/2 = 0.14). 
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Section 3.4 uses these, and subsequent eigenvectors, in spectral expansions of 
the dynamics within the ith patch in order to discover the effects of different 
mesoscale coupling procedures (Section 4). 


3.3 Action regions and core average over part of a 
patch 


Let’s now consider the case where the core and action regions have half-width 
a > 0. As in Section 3.2, we evaluate numerical examples of eigenvalue 
and eigenvectors using Matlab and from these we determine the analytic 
forms, which we confirm by substitution into the matrix equations (9). The 
eigenvalues of matrix equations (9) are 


Ak = 



cos[7rl lc /2(n — a)]}, k = 0,1,..., 2(n — a — 1), 
cos[7tl k /(2a + 1)]}, k = 2(n— a) — 1,... ,2(n— 1). 


The k < 2(n — a — 1) eigenvalues are similar to those for the a = 0 case, 
except that the half-width of the patch is now effectively (n— a) rather than n, 
thus reducing the buffering of the macroscale solution from ntou-a. For 
k = 0,1,..., 2(n — a — 1) the wavenumbers are the same as the a = 0 case, 


lie = 


k + 1 + (—l) k/2 (2£/7t— 1) 

k+1 


for even k, 
for odd k. 


(15) 


For k = m + 2(n — a — 1) where m = 1,2,..., 2a, the wave numbers are 

Vn+2(n—a—1) = 2 [m/2]. (16) 

The last 2a eigenvalues are a equal pairs. If we set a = 0, then the above 
eigenvalues reduce to those given in equation (10). Figure 7 plots scaled 
wavenumbers against eigenvalues for patch half-width n = 20 and action and 
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Figure 7: Microscale eigenvalues plotted against wavenumbers which are 
scaled by the reduced width of the patch, l k /2(n — a), for k < 2(n — a — 1) 
(blue) and scale by the width of the core l k /(2a + 1), for k > 2(tl — a — 1 ) 
(red). The patch half-width is n = 20, the action and core half-width is 
a = 5 and cost = 0.91 . 
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core half-width a = 5. For this case there are 2(n — a — 1) + 1 =29 unique 
eigenvalues and a pairs of equal eigenvalues. 

The right eigenvectors for all sub-patch coordinates j and k = 0, 1 ,..., 2(n — 
a — 1) are 


(Vk)j 


'll — a) sin i\ 1/2 cos[jl k 7r/2(u 
(n — a )^ /2 sin[jl k 7t/2(n — a)] 


a)] for even k, 
for odd k. 


(17) 


If we set a = 0 , then the above eigenvectors reduce to the right eigenvectors 
given in equations (12) and (13). For k = m + 2(n — a — 1) where m = 
1,2,..., 2a, the right eigenvectors for all sub-patch coordinates j are 

_ f (2a + 1 )^ 1/2 cos[jl k 7t/(2a + 1 )] for even k, 
k ’ [ (2a + 1 )^ 1/2 sin[jl k 7t/(2a + 1)] for odd k. 


Despite there being a pairs of equal eigenvalues for k > 2(n — a — 1 ), the 
associated right eigenvectors are linearly independent since for even k the 
eigenvectors are even in j and for odd k the eigenvectors are odd in j. 2 

The left eigenvectors are considerably more complex than the right eigenvec¬ 
tors. To define the left eigenvectors we first define some new functions. For 
k = 0,1,... ,2(n — a — 1), 


(wic)i = (-1) rtk 1)/21 {2sin[(2a+ 1)l k 7i/4(n- a)]} 


-l 


x 


cos 


l k 7T 


4(n — a) 


— cos 


(2j + 1)l k 7t 
4(n — a) 


2 It is possible, for some particular n and a, for there to be odd ki < 2(n — a — 1 ) and 
odd k 2 > 2(u — a — 1) such that A k| = A k2 which implies v k| = v k2 . For example, such 
a case occurs for n = 4 and a = 1 when k] = 3 and k 2 = 5 . Such repeated eigenvectors 
mean that we do not have a complete basis. To form a complete basis we should replace 
repeated eigenvectors with generalised eigenvectors [4, e.g.], but we avoid this complication 
by avoiding those n and a which result in repeated eigenvectors. 
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[(n — a) sin £] ^ 2 for even k, 
(n — a)~ 1/2 for odd k, 


(19) 


and for k = m + 2(n — a — 1) where m = 1,2,..., 2a, 


( Wk)j = (- 1) r(m “ 2)/21 


cos 


l k 7T 


2(2a + 1' 


— cos 


( (2j + 1)lk7t\ 
\ 2(2a + 1) ) 


(2a + 1 ) 1/2 {sin[(n — a)l k 7t/(2a+ 1)]} \ 

(2a + 1 ) 1/2 {cos[(n — a)l k 7i/(2a+ 1)] — cosf}^ 1 , 


( 20 ) 

for odd m, 
for even m. 


We also define, for k = 0,1,..., 2(n — a — 1) and |j| = 0,1,..., n — a, 



[(n — a)sin£] 1/2 sin[f — (—1 ) k/2 |j|l k 7x/2(n — a)], even k, 
(n — a)~ 1/2 sin[jl k 7r/2(n — a)], odd k, 


and for all other possible values of k and j we set (z k )j = 0. 

The left eigenvectors for sub-patch coordinate |j| < n — 1 are, for even k, 


(^k)j — (z k )j + < 


' -2 cosf(w k ) a _|j|, 

0, 

( w k)jj|—n+2a i 
, (W k ) n _|j| , 


o < lj| < a, 
b < |j| < ri —2a, 
n — 2a < |j| < n — a, 
rt — |j| <ru — 1 , 


and for odd k, 


( 22 ) 


[0, 

(z k )j = (z k )j + < sgn(j)(w k ) B |_ n+2a , 
[sgn(j)(w k ) n _|j|, 


0 < |j| < n — 2a, 
n — 2a < |j| <n—a, 
n—a<|j|<n—1. 


(23) 


In all cases, (z k )± n = (z k )-t-( n _i) . These z k reduce to the left eigenvectors 
given in equations (12) and (13) for the special case a = 0. 
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Equation (19) is undefined when k = ki < 2(n — a — 1) is odd and (2a + 
1)l kl 7t/4(n — a) = S7T for integer s, or equivalently, when l kl 7t/2(n — a) = 
l k2 7r/(2a +1) for k 2 = m + 2(n — a — 1) and m = 2s — 1 . From the 
limits on ki it can be shown that s = 1,2,..., a so m = 1,3,..., (2a — 1). 
Therefore, equation (19) is undefined if for odd k = ki < 2(n— a — 1) there 
exists odd k 2 > 2(n — a — 1) such that A kl = A k2 ; a scenario which we 
identify as a degenerate case which requires generalised eigenvalues and is 
avoided. Similarly, equation (20) is undefined when k = k 2 > 2 (tl — a — 1) 
is odd and (n — a)l k2 7t/(2a+ 1) = S7t for integer s, or equivalently, when 
l k2 7t/(2a + 1) = l kl 7t/2(n — a) for ki = 2s — 1 . From the limits on k 2 it can 
be shown that s = 1,2,..., (n — a — 1) so ki = 1,3,..., 2(n — a — 1) — 1 . 
Therefore, equation (20) is undefined if for odd k = k 2 > 2(n — a — 1) there 
exists odd ki < 2(n— a— 1) such that A k2 = A kl , which is, again, a degenerate 
case to be avoided. 

Figure 8 plots right and left eigenvectors associated with the four lowest 
magnitude eigenvalues with k < 2(n — a — 1), for patch half-width n = 20 
and a = 5, and Figure 9 plots right and left eigenvectors associated with the 
four lowest magnitude eigenvalues with k > 2 (tl— a — 1), with the same patch 
geometry. These plots are typical of any n and a ^ 0, provided we have a 
complete basis of eigenvectors across sub-patch coordinates |j| = 0,1,..., n— 1 
(so do not require generalised eigenvectors). The shape of the right eigenvectors 
are not remarkably different from those for a = 0 shown in Figure 6, with the 
main point of difference being the frequency of the sinusoidal eigenvectors 
which are increased by the effective reduction of patch size from n to n — a. 
In contrast, the left eigenvectors for a > 0 are significantly different to those 
for a = 0. For nonzero action regions and core, z k with odd k < 2(n— a — 1) 
appears smooth about j = 0, unlike z k for zero action regions and core. For 
k > 2(n — a — 1) the left eigenvectors with odd k are only nonzero inside the 
action regions whereas for even k the left eigenvectors are only nonzero inside 
the action regions and the core. 
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Figure 8: Right and left microscale eigenvectors associated with the four 
lowest magnitude eigenvalues with k < 2(n — a — 1) (k = 0,1,2,3 for blue, 
green, red and cyan, respectively) with the same parameters as Figure 7, 
patch half-width n = 20, action and core half-width a = 5 and cost = 0.91 . 


Thursday 9 th April, 2015 






3 Microscale sub-patch dynamics 


27 




Figure 9: Right and left microscale eigenvectors associated with the four 
lowest magnitude eigenvalues with k > 2(n — a — 1) (k = m + 2(n — a — 1) 
with m = 1,2,3,4 for blue, green, red and cyan, respectively) with the same 
parameters as Figure 7, patch half-width n = 20, action and core half-width 
a = 5 and cos£ = 0.91 . The right eigenvectors are only plotted over the 
patch core, |j| < a, as this shows at least one complete period. The left- 
eigenvectors are plotted over the entire patch, |j| < n. 
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3.4 Microscale field for continuous coupling 


Given the homogeneous microscale solutions of the previous subsection, we now 
explore a spectral representation of the solution within a patch with forcing 
by general coupling with neighbouring patches. The microscale field solution 
of ODE (7) within the ith patch, which has continuous and instantaneous 
coupling with nearby patches defined by coupling conditions (3), is 


pt 


2(n—1) 2(n—1) 

u(t) = Y_ e Akt v k z[Bu(0) + fit) + Y_ v k zl f(t V k(t-t,) dt' 

k=0 k=0 Jo 

= T(t)Bu 0 + fit) + T(t) * fit), (24) 


where u = (u_ n+ iN,..., u u+ in) describes all microscale fields Uj+iN with 
|j| < n within the ith patch, and the convolution in the last term is defined 


as 


2(n—i) t 

T(t) * fit) = Y_ T k fit V 1 ^'’ dt', 


k=0 


(25) 


with (2rt + 1) x (2n + 1) matrices T k = v k z^ and state transition matrix 


2(n—1) 

T(t)= Y_ eAktT K- (26) 

k=0 

This solution is confirmed via direct substitution into ODE (7) and with the 
identity, proved in Appendix A, that 

2(n—1) 2(n—1) 

T(0) = Y_ Tk = Y- ^ = B + A, (27) 

k=0 k=0 

where all elements of the (2n + 1) x (2rt + 1) matrix A are zero, except 


Aj,± n = £j >±n and A ±njj = £ ±nij for |j | = 0,1,..., u. (28) 
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Thus, on substituting solution (24) into ODE (7) we use 

2(n—1) 

B Y_ Tkf(t) = BAf(t) = (0, Vn, 0,..., 0, ft, n , 0) = f(t) + £f(t). (29) 

k=0 

Identity (27) is also useful in confirming solution (24) has the correct initial 
condition. 

The general forced solution (24) acts as the reference for reporting errors due 
to infrequent mesoscale coupling. 


4 Infrequent mesotime coupling 

This section approximates the exact analytic solution (24) of ODE (1) within 
one patch with coupling conditions (3) by replacing the continuous time 
coupling vector f(t) with a mesoscale coupling vector. With mesoscale 
coupling we only evaluate the coupling vector f(t) at fixed time intervals 
t = m5t where the meso-time interval 5t is much larger than the time-step 
of the microscale computation. We assume we know the first Q terms of 
the Taylor series in time of the coupling f (t) exactly at times t = m6t for 
m = 0,l,...,M, but at no other points in time. If we only know the values 
of f, and not their time derivatives, then parameter Q = 1 ; but our analysis 
applies for general Q. Further, we assume the computation in a patch up to 
time t = m6t only depends on f (t 7 ) evaluated at t' < m6t . Thus we know 
the (Q — 1 )th order Taylor series of f(m5t +1) about t = 0 for 0 < t < St 
and we apply coupling conditions (6). We first consider the initial time step 
with m = 0 where the coupling f(t) and its derivatives are only known at 
t = 0 . By homogeneity in time, the analysis extends straightforwardly to a 
general number of meso-time steps. 
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4.1 Initial time step 


Consider the solution (24) at t = St multiplied by matrix B, 


2(n—1) 2(n—I) 

Bu(6t) = B [a. k T k Bu 0 + BT k J k , 

k=0 k=0 


where q k = e Ak5t and 



f(t / )e Ak(6t ~ t/) dt'. 
o 


(30) 


(31) 


The multiplication by B in the above solution removes the fields at the 
endpoints of the patch, u± n , which can always be evaluated from the coupling 
conditions (3), (5) or (6), if required, provided all other held components 
are known. The aim is to approximate J k in (30) using the (Q — 1)th order 
Taylor series expansion of f(t'), as used in coupling conditions (6), and also 
to provide the error of J k due to this Taylor series approximation. 

Integrate (31) by parts, 



T“fi0)(p k — 1) + T~ 

Ak A k 


o 


1)dt', 


(32) 


where a constant of integration is chosen so that only f (0), not f (5t), appears 
in the integrated part. The first term in equation (32) is exactly what is 
obtained by substituting the zeroth order Taylor series of f(t') about t' = 0 
into integral (31) (i.e., replace f(t') with f(0)) and thus the integral part 
of (32) is the error of the zeroth order Taylor series approximation. Then, 


2(n—1) 


Bu(6t) = Y_ BT i< M-kBuo + f(0)(p k — 1 )/A k 


+ BR, 


(33) 


k=0 
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with remainder vector R = (R_ n , R_ n+ i,..., R n ) where, for \j \ — 0,1,..., (n — 

D, 


2(n—1) 


«.= L 


i 


k=0 

2(n-1) 

L 

k=0 


■6t 


BT k f(t')l (e Ak(6M ' ) -1)dt' 


(Vk)j 


p6t 


0 L 


(Z k )_ n f-u(t') + (z k ) n f n (t')l - 1)dt' . (34) 


Since the microscale fields in the action regions are constrained by the cou¬ 
pling conditions (5), the remainders at j = ±rt are dependent on the other 
remainders in the action regions, 


±(n-1) 

R±n = — Rj • 

j=±n—2a 

In the mesoscale coupling scheme with coupling conditions (5) we use 


(35) 


2(n—1) 

Bu(6t) - Y_ BTk 

k=0 


p k Bu 0 + f(0)(p k — 1 )/A k 


(36) 


with the error known to be precisely BR. 

The solution (32) for J k is for one integration by parts which is equivalent to 
approximating f(t') by a zeroth order Taylor series (Q = 1), as in coupling 
conditions (5). Generalising to Q integrations by parts, 


q—i 


h=Y_ fVl «» V bk - X 6t? A£/p! 


q=1 

1 

+ 4? 


p=0 


p6t 


f^t') 


Q-1 


Lll! A P( t ' _ 5t)P 


p=0 


P- 


dt', 


(37) 
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where we choose constants of integration so that only f(0) and its derivatives, 
not f(5t) and its derivatives, appear in the integrated part. In this case 


2(n—1) 

Bu(6t) = Y_ BT >< 


k=0 


p k Buo+ y~ 


q=i 


f»(0) 




+ BR, (38) 


with components of the remainder vector 


2(n-1) 


1 


R <= L a qj 

k=0 A k 


St 


[BTkf^ltOJi 


b 1 r_i’ip 

,A k (6t-t') _ y 1- 5t) p 

Z— D 

p=0 K 


dt' 


h 1 rSt 00 r -I'm 

Y_ 7q (^[(&U^) + (4)nf2(tO]i;-^A£(t'- 6t) p dt' 

Av Jo p= q V- 


k=0 /4 k 
_ 2(n—1) 


L 

p=0 


L 


AE 


L 

p=0 


k=0 
2(n—1 


(Q + p)! 


L 

k=0 


A? 


(Q+P)! 


(v k )j(z k )_n 


(vk)j(4) T 


p6t 


fS n (t')(5t —t') Q+p dt' 


p6t 


f?(t')(6t — t')^^ dt', 


(39) 


for |j| = 0,1,..., (n — 1), and where R ±n satisfy (35). 


On substituting the (Q — 1)th order Taylor series of f(t') into integral (31), 
a general term in the expansion, for q = 0,..., (Q — 1), is 


1 

q! 


rSt 


dt' = f q (0)A k (q+1) ( Hk - X. 5tPA k/P ! ) > (40) 

p=0 


which describes all terms in the first line of the expansion of J k in (37). 
Thus, this first line of (37) is the approximation of J k with mesoscale coupling 
conditions (6) obtained from the (Q—1 )th order Taylor series expansion of f(t) 
and the second line of (37) is the error of J k due to the Taylor series expansion. 
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In the mesoscale coupling scheme, when we know (Q — 1) >0 derivatives 
of f(t) at t = 0 we improve on the mesoscale coupling conditions (5) and 
approximate field solution (36) by using coupling conditions (6) which provides 
the approximate field solutions 


2(n—1) 

Bu(St) ~ Y_ 


BTv 


k=0 


p k Bu 0 + ^f Vl (0)A k 


q-1 


Bk 


£ StfApp! 


q=l 


p=0 


(41) 


with the error known to be precisely BR. 


4.2 Multiple mesotime steps 


After M mesoscale time steps of size St, solution (24) gives 

2(n—1) 2(n—1) M—1 

Bu(MSt) = B Y_ Bi^T k Bu(0)+ BT k ^<— ] J km (42) 

k=0 k=0 m=0 


where the integral from time zero to M5t is converted into M integrals from 
zero to St, 

p6t 



f(t'+ m6t)e Xl|SM '> dt'. 


(43) 


0 


The integral J km is a generalised version of the integral J k in equation (31). 
After Q integrations by parts, J km is similar to J k in equation (37), but with 
f(0) and f(t') replaced with f(mSt) and f(t 7 + mSt), respectively. Thus we 
obtain 


2(n—1} 

Bu(MSt) = Y_ BT >< 

k=0 


Q 

p^Bu(O) + Y_ 

q=1 





q-l 

Bk- Y_ 

p=0 



+ BR 

(44) 
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where 


M-1 

m=0 


the components of the remainder vector for |j | = 0,1,..., (n 


1) are 


(45) 


2(n—1) 




r§ t 


A Q 
k=0 A k 




d 1 (■ _I 

,A k (6t-t') _ y 1—L A P( t ' - 8t) v 

^— n! 

p=o p 


dt', (46) 


and R± n satisfy (35). 

The solution over M mesoscale time steps (44) is the same as the solution over 
one mesoscale time step (38), but with derivatives of the coupling vector f(0) 
replaced by derivatives of the M time steps coupling vector f M (0). The 
only difference between the remainder vector R over M mesoscale time steps 
and over one mesoscale times step, is that for M, equation (46) integrates 
over f^t') and for one time step, equation (39) integrates over . In 

the following error analysis we only consider a single mesoscale time step, but 
these results are readily adaptable to multiple mesoscale time steps. 


5 Error analysis 


The remainder vector R (39), or (46) for multiple mesoscale time steps, tells 
us the extent to which coupling errors penetrate into the core of the ith patch 
from the action regions. Here we are only concerned with errors arising 
from evaluating the patch coupling f(t) at mesoscale time intervals in the 
mesoscale coupling conditions (5) or (6), not with errors which are due to, 
for example, the interpolation between patches. Ideally, the remainder Rj 
should be minimal within the patch core (i.e., for |j| < a) since it is the 
microscale fields iq within the patch core which determine the macroscale 
held Ui(t). Errors near the patch edges (i.e., j ~ ±n) need not be small. 
Section 5.1 considers how Rj varies across the patch, |j| = 0, and 
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with varying the number of terms Q in the Taylor series expansion in the 
coupling conditions (6). Section 5.2 then looks at the upper bound of the 
error of the macroscale held Uj.(St) due to mesoscale coupling. This analysis 
considers one mesoscale time step as the base to predict multiple steps. 


5.1 Error penetration in patch 

Defining fO(St) = max t ' e [o,st] [|fi]L n (t , )|, l"fpn(l / )l] and rearranging equation (39) 
gives the scaled upper error bound for the remainder when |j| < (n — 1), 


IRi 


< 


st Q+1 


2(n—1) 


fQ(5t) ~ (Q+l)! "I” 1 E 1 ( 1 ; Q + 2;A k 5t)| — Rj max , 

(47) 

where the sum over p is rewritten as a confluent hypergeometric function [13] 


L 


a p r 6t 5tQ +1 

k (6t - t') Q+T, dt' = — ; tF, (1; Q + 2; A k St). (48) 


(Q + P)! 


(Q + l 


For j = ±n, from equation (35), 


±(n-11 


k ". max ^ k j ma x 

j=±n—2 a 


(49) 


Figure 10 plots the upper bound of the remainder Rj m ax f° r a range of Q 
and St = 0.5 . Since the remainder is symmetric about j = 0 , only the 
j > 0 components are shown. Figure 10 shows that the remainder is very 
small in the core of the patch, as required. For St ~ 0.5 the upper bound is 
approximated by 

R jma x - (2St) Q+1 10- Q - (1+0 - 025Q/St)(u - 1 -’ ) (50) 

until numerical error dominates around Rjmax ~ 10~ 13 ~d. With constant 
patch half-width n and |j| < n, the core half-width a has little affect on the 
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remainder Rj max and, until numerical error is dominant, there is nothing to 
distinguish the remainders with different a, provided the mesoscale time St 
and number of terms in the Taylor series expansion Q are fixed. At the 
patch edge |j| = n the remainder upper bounds for core half-widths a > 
0 are distinctly different from the a = 0 case. For a > 0, in (49) the 
R±(n—i) max term dominates the sum and R ±nmax - R±( n -i)max • For a = 0 the 
sum in equation (49) vanishes and thus R± nm ax = 0 . 

The inset of Figure 10 shows that varying cost, which is a function of the 
ratio r = (n — a)/N and the coupling strength y, does not affect the upper 
bound of the remainder. Furthermore, as implied by the approximation (50), 
the upper bound of the remainder is dependent on the distance from the 
patch edge, that is n — j , rather than the patch half-width n alone, and Rj max 
decreases exponentially with n — j. Therefore, larger n does not in general 
produce smaller Rj max ; however, it does produce wider patches which have 
smaller Rj max within the centre core of these patches. In macroscale modelling 
the most important error measurement is the error in the macroscale held, 
and the macroscale held is obtained by averaging over microscale helds in 
the patch core. Therefore, for the system considered here the errors in the 
macroscale held will generally be minimised for larger patch half-widths n 
and smaller core half-widths a. 

We expect similar results for other microscale systems in the same universality 
class as (1). 

5.2 Macroscale solution 

The macroscale held value obtained from the ith patch, Ut (St), is the average 
of the microscale held Uj + tN(5t) in the patch core, as shown in equation (2). 
The error of this macroscale value is the average of the remainder vector 
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Figure 10: The solid lines describe the upper bound of components of the 
remainder vector Rj max in a patch with patch half-width n = 20 , cos l = 0.91 , 
mesoscale time 6t = 0.5 and Q = 1,3,5, 7. Lines with the same colour 
have the same Q but all possible core half-widths a = 0,..., 19; these 
different core half-widths are only distinguishable in the centre of the patch 
j < 10, and on the very edge of the patch j = 20 for the a = 0 case. The 
dashed lines approximate the coloured curves with a simple power function. 
Numerical error dominates at low Rj max . Inset: the remainder vector Rj max 
on the same scale as the main plot with n = 20, 6t = 0.5, Q = 1 , and 
cos i — 0.65,0.75, 0.85,0.95 . Lines with the same colour have the same cos i 
but all possible core half-widths a = 0,..., 19 and are only distinguishable 
below Rj max - 10~ 13 where numerical error dominates. 
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components Rj in the patch core: 


E = 


1 


2a + 1 


L*. 

j=—a 


1 


2a+1 


L 

p=0 


2(n—1) 




L (QTrtiLw.a) 


k—0,even 


j=-a 


rSt 


X 


[fl n (t') + fS(t')](6t — t') Q+p dt' 


(51) 


where only the even eigenvalues and eigenvectors contribute. 

Similarly to the scaled upper bound of the remainder in (47), the scaled upper 
bound of the error is 


(2 a _i_ 1 'ip k+Q+i 2(n_1) a 

f Q (S t ) - (Q + 1)! H T. l(vk)j(zk)niFi(1;Q + 2;A k 5t)| = E max . 

1 ^ ’ k=0,even j=—a 

(52) 

In this error bound we scale by (2a + 1) because fQ is of the order (2a + 1) 
due to scaling of the coupling conditions (6) and, for example, (4). 

Figure 11 plots the scaled error bound E max for Q = 1 (communicating only 
function values f± n (0) and no derivatives as in coupling conditions (5)) and 
a range of mesoscale time-steps 6t. This figure shows that the upper bound 
of error E max is a function of reduced patch half-width n — a, rather than 
individual values of n and a, and decreases as n — a increases. So, the error 
is smaller for large n and small a, as predicted at the end of Section 5.1. 
Section 5.1 also showed that the remainder upper bound Rj max is independent 
of cos £, and thus independent of ratio r = (n — a)/N , and we conclude that 
the error upper bound E max should then also be independent of cost and r. 
This is a surprising result as it implies that the error is independent of the 
ratio between the microscale and macroscale lattice spacings N = H/h.. We 
interpret this to mean that it is most important to make the reduced patch 
half-width n — a (or buffer width) only large enough to capture a significant 
portion of the microscale dynamics, and then macroscale modelling may 
proceed to any scale. The interpretation of ‘significant portion’ will depend 
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on the nature of the microscale model. Of course, here we are exploring 
an error bound scaled by fQ(6t); some r dependence will be present in the 
absolute error when, for example, we use of the form (4). 

Even for only one term in the Taylor series expansion of the coupling condi¬ 
tion (5), Q = 1 , the upper bound of error E max plotted in Figure 11 decreases 
significantly with increasing n — a or decreasing 6t. Based on Figure 10 we 
expect larger values of Q to produce smaller errors for macroscale field Ut(t), 
but in a simulation, increasing Q might not be practical. This is because a 
given Q requires the first Q — 1 temporal derivatives of the coupling vector f, 
and this information may not be practically available. Whether it is more 
practical to increase n — a or decrease 5t in order to reduce errors is dependent 
on the problem being considered and the architecture of the parallel com¬ 
puter running the simulation. For example, in cases where there is periodic 
microscale detail, choosing n — a to be multiples of the period gives more 
accuracy [6] and so it might be best to choose an optimal n — a and then con¬ 
sider reducing 5t. In general, while increasing n — a will increase processing 
times, decreasing 6t will increase the amount of data communication. In the 
case of a supercomputer with a large number of processors, maximising n — a 
will take advantage of the processing power while choosing 5t large enough 
will avoid limitations associated with interprocessor communication. 


6 Two dimensional Ginzburg-Landau 
numerical simulation 

To demonstrate that mesoscale coupling is effective in dimensions higher 
than one and in more complicated problems than the simple diffusion of 
ODE (1), we use patch dynamics with mesoscale coupling to simulate the 
complex Ginzburg-Landau equation [1] on a discrete two dimensional square 
microscale spatial lattice with lattice spacing h: 

<,j y W = (1 +^)[u jx+1)jy (t) +U jx _ 1)jy (t) -2uj x> j a (t)] 
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Figure 11: The upper bound of the error E max for Q = 1 and cost = 0.91 , 
a range of mesoscale time steps St and over several reduced patch half¬ 
widths n — a . The calculated range of the patch and core half-widths was 
for all 4 < n < 20 and 0 < a < n, respectively, but since the error E max 
is extremely small for n — a > 11 these cases are not plotted. Curves with 
different n and a are identical when they have the same reduced patch half¬ 
width n — a and mesoscale time step St (and are thus plotted in the same 
colour), until numerical error dominates the calculation below E max ~ 10~ n . 
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+ (1 + ia) [Uj x , jy +i (t) + Uj X)jy _i (t) - 2u jx , jy (t)] 

+ u ix,j y W-n + i PHx,jJ t )K,j y (t)l 2 , ( 53 ) 

for constant, real a and (3. The complex Ginzburg-Landau equation contains 
two dimensional diffusion terms with an additional nonlinear part. We set 
a = 1 and (3=2; this choice of parameters should, after some time, produce 
plane wave solutions. The initial conditions specified were sinusoidal with 
amplitude 0.5 plus a real, normally distributed random number with mean 0 
and standard deviation 0.8 at each microscale lattice point. 

The construction of a two dimensional macroscale lattice with lattice spacing H 
and two dimensional patches is very similar to the one dimensional case 
discussed in Section 2. As in the one dimensional case, we assume that 
N = H/H is an integer. We construct the square macroscale lattice with 
general lattice point (X ix , Yj. ) = (i x H,iyH) for integers t X) y . Centred about 
each macroscale lattice point (Xi x , Yi y ) we construct the (i x , i y ) square patch of 
width 2nh which does not touch or overlap neighbouring patches. Again, like 
the one dimensional case, integer n is the patch half-width and to ensure the 
patches do not overlap n/N < 1 /2 . The microscale fields on the (ix, iy) patch 
are ib x +i x N,j y +iyN for microscale sub-patch index j Xj y = 0, ±1 ,..., ±n. In the 
patch dynamics numerical simulation the Ginzburg-Landau equation (53) is 
solved for all microscale values within a patch, excluding those on the patch 
edges, j x , y = 0, ±1,..., ±(n — 1). 

We could define the macroscale fields in terms of some average over microscale 
fields in the centre or core of the patch (i.e., average over fields in a square core 
with core half-width a), as we did in the one dimensional case (2). Similarly, 
we could define patch coupling conditions in terms of some average over 
microscale values on the patch edge, as in equations (3), (5) or (6). However, 
Section 5 showed that mesoscale errors are a function of the reduced patch 
half-width n — a rather than patch half-width n so here we simply choose the 
two dimensional analogue of a = 0 where we act and sample the microscale 
field at patch edge points and patch mid-points to define patch coupling 
conditions and macroscale holds, respectively. Thus, the macroscale held is 
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the microscale field in the centre of the square patch [31], 

(t) = ttixN^yN (t) , (54) 

and the patch coupling conditions with mesoscale coupling St constrain each 
microscale field on four patch edges, 


u± n +i x N,j y (mSt + t) = U ix , iy (mSt + t) cos£ ±n , jy + f(i x ,i y ),(±n,j y )(mSt) , 
Uj X)±n +i,N(rnSt + t) = u ix)iy (mSt + t) cos £ jx , ±n + f ( i x ,i y ), ( j x ,± n ) (m6t), (55) 


for nonnegative integer m and j x>y = 0, ±1,..., ±(n — 1). These coupling 
conditions are analogous to the Q = 1 one dimensional coupling conditions (5). 
Both cos£j X) j y and f[i Xl i y ),(j Xl j u ) are functions of the patch coupling strength y 
and the ratio r = nfi/H and describe interpolation from the centre of the 
(i x ,iy) patch and across the un-simulated space between the (i x ,iy) patch 
and adjacent patches [31]. Coupling conditions are not required for microscale 
fields on the patch corners where both j x = ±rt and j y = ±n since these 
microscale fields are not required when solving the Ginzburg-Landau equa¬ 
tion (53) for the microscale fields within the patches but not on the patch 
edges. 

For the numerical simulations presented here we use nearest neighbour cou¬ 
pling derived from Taylor series expansions about the microscale points on 
the patch edge [31]: 


cos£ jx,j y 

%x,iy),(jx,jy) W 


1 -( r fx- r t.)T, 


hr 


i r ixy[( r ix ± 1)U ix+1 ,i y (t) + (r jx =F l)U ix _i )iy (t)] 

+ Ky[(r jy ± 1 ) U ix,i y +i (t) + (r jy T 1 )U ix , iy ^ (t)], (56) 


for tj x = j x h/H and p y = j y h/H where j Xjy = 0, ±1,..., ±n and r ±n = ±r. 
We define a square periodic domain of width 20 with microscale lattice spacing 
h = 0.25. We set macroscale lattice spacing H = 5, which fits 4x4 patches 
across the domain and gives N = H/h = 20. We choose a patch half-width 
n = 6 so that r = 0.3 and choose mesoscale coupling 6t = 0.2. 
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Figure 2 compares the real components of the microscale fields within patches 
for continuous time coupling and with mesoscale coupling at two times, 
t = 0.04 and t = 0.4 ; there is little difference between these two cases. 
The large random component in the initial condition decays rapidly and is 
substantially reduced even at small time t = 0.04 . By t = 0.4 the simulation 
is smooth. Figure 12 compares the real and imaginary part of the macroscale 
fields U ixjiy obtained with continuous time coupling and mesoscale couplings 
St = 0.2,0.1 . The St = 0.2 case produces a reasonable result, but some 
macroscale fields are slightly inaccurate. Better results are obtained when 
the mesoscale coupling time-step is reduced to St = 0.1 . 


7 Conclusion 

We analysed a patch dynamics macroscale modelling scheme which is adapted 
for massive parallelisation and exascale computing by limiting the transfer 
of data between processors, assuming that one processor only calculates the 
dynamics of one patch. Rather than transferring information concerning 
coupling conditions at every microscale time step, we propose limiting the 
data transfer to mesoscale time-steps St. This method does not require any 
preprocessing and makes minimal use of stored data. Limiting the transfer and 
size of data addresses several major hurdles facing the development of exascale 
computing, specifically, slow data transfer speeds compared to processor 
speeds, the energy cost of data transfer and memory size limitations [36, 37,11]. 
If one processor was to evaluate the dynamics of a small number of adjacent 
patches, then coupling data should be updated at microscale times for the 
patches evaluated on that one processor, but coupling data transfers at 
mesoscale times should be maintained for patches evaluated on different 
processors. 

Section 5 found that errors arising from mesoscale coupling are controlled 
by the reduced patch half-width rt — a and the mesoscale coupling time St, 
with larger n — a and smaller St producing smaller errors. However, these 
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time t 


Figure 12: Real and imaginary parts of macroscale fields U lxi i H from all 
16 patches with patch half-width n = 6. Compare continuous time coupling 
(solid lines) and mesoscale coupling plotted at the mesoscale time steps 
6t = 0.2 (open circles) and 6t = 0.1 (dots). 
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error predictions only compare patch dynamics with continuous time coupling 
to patch dynamics with mesoscale coupling; they do not consider errors 
inherent in the patch dynamics scheme due extrapolation across large regions 
of un-simulated space, and these additional errors are also functions of n — a. 
Previous work analysed the error of patch dynamics with continuous time 
coupling relative to the known exact solution and found that larger n — a 
are not always better when the model has rough microscale detail [6]. If 
the microscale detail has some periodicity, then optimal solutions require 
that rt — a is chosen such that the periodicity exactly divides n — a. Thus, 
adjusting the reduced patch half-width rt — a for improvements in mesoscale 
coupling patch dynamics is not practical if rt — a is already constrained 
by the symmetry of the microscale model. In these cases, to reduce errors 
for mesoscale coupling it is advisable to fix ri — a at the optimal solution 
determined from the symmetry of the microscale model and only reduce 5t. 

Section 1 briefly discusses the need for fault management and resiliency in 
exascale computing [36, 37, 11]. While the proposed patch dynamics scheme 
makes no allowances for failure, we expect that some modifications to the 
scheme will address this issue. Future work may develop coupling conditions 
for patch dynamics which are dependent on variable mesoscale times steps 5t, 
so that some delay in data transfer is accounted for in the algorithm. In 
addition, this variable time step may be permitted to extend to infinity, thus 
accounting for cases where the data never arrives, due to, say, complete failure 
of a processor. 

The mathematical analysis presented here only considers a simple one dimen¬ 
sional microscale diffusion model, but the scheme is readily modifiable to 
patches in two or more spatial dimensions and more complex models [31]. The 
numerical simulation in Section 6 show that mesoscale temporal coupling is 
effective for the complex two dimensional nonlinear Ginzburg-Landau model, 
and reducing the mesoscale coupling time 5t produces more accurate results, 
in agreement with the one dimensional analysis. More work need to be done 
to develop the full implementation of patch dynamics for mesoscale coupling, 
that is, with patches in both space and time [19] illustrated in Figure 1. 
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A Proof of identity (27) 

Here we prove the identity (27) which defines the state transition matrix T(t) 
at time t = 0. The proof holds for C defined in Section 3.1, provided the 
matrix equations (9) produce k = 0,1,..., 2(rt— 1) linearly independent right 
and left eigenvectors, v k and z k , which satisfy the normalisation condition 
v k Bz= 5 kk / for all k, k' = 0,1,..., 2(n — 1). The set of all v k or all z k form 
a complete basis which spans the subspace of the microscale held vector on 
one patch u satisfying the coupling conditions (3), (5) or (6). 

Define some vector in the subspace of the microscale held vector u in terms 
of the right eigenvectors 

2(n—i) 

s = B Y_ C DT, (57) 

k=0 

for arbitrary coefficients c k . The hrst and last components of this vector are 
s± n = 0, but, since the right eigenvectors form a complete basis for sub-patch 
coordinates j = —n + 1, —n + 2,..., n — 1, all other components of s are 
completely general. Using the normalisation condition of the eigenvectors, 

2(n—1) 2(n—1) 

BT(0)s = B Y_ v k z k Bc k /V k / = B Y_ c Wk = s, (58) 

k,k'=0 k=0 

and since the form of s is arbitrary for all but its hrst and last components 
we conclude that BT(0) = B + M where the matrix M is such that Ms = 0. 
Since s is a general vector, except for s± n = 0, the only nonzero elements 
of M are in its hrst and last columns. In addition, since all elements in the 
hrst and last rows of BT(0) must be zero, we also know that all elements in 
the hrst and last rows of M must be zero. We define some matrix A which 
satishes M = BA and, given the form of M, conclude that the only nonzero 
elements of A are in the hrst and last rows and the hrst and last columns. 
Therefore, 

T(0) = B + A . (59) 
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We now show that the matrix A is of the form given in equation (28). 

From equation (9) we obtain [z^C)± n = 0. Thus, for ) = —n + 1,—n + 
2 ,..., n — 1 , 

2(n—1) 2(n—1) 

rr(0)£]j,±n = o = Y- [wzk£] j)±n = Y_ (vk)j(Zk£) ±n 

k=0 k=0 

= T Aj,— n£— n,±n T • (60) 

Therefore, for j yf ±n, using the definition of C in Section 3.1, 

^•j,±n — Sj,=t(n—1) — ^j,±n • (61) 

Similarly, from equation (9) we obtain (£vk)± n = 0, so, for, j — — n, — n + 

2(n—1) 2(n—1} 

[£T(0)]± n j = 0 = Y_ [^V k Zk]±n,j = Y (^W)±n(Zk)j 

k=0 k=0 

= £±n,j + _ n A_n,j + £-j- n)n A n) j . (62) 

For j yf ±n and using equation (8), 

A±n,j = £± n ,j + ~ ^±n,j • (63) 

Finally, since [£T(0)] ±rvn , = [£T(0)]± n _ n = 0, 

n 

^ Tlin^A^n = £±n,n A n>n T 22±rL,(rL—1) —1) , n = 0 , 

l=—n 
n 

^ '£±n,l'A-l,—n f-'±.n,— n^-—■n,—'n T 21±n,—(n— 1 )A_(n—1),—■n 0 , (64) 

l=—n 

so that 

^■±n,±n = 1 = l-'±n,=tn and A^^ip^ = 0 = 21±n,=pn • (65) 

Therefore, with equations (61), (63) and (65) we have shown that the matrix A 
is of the form given in equation (28). 
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